The Growing Correlation Length in Glasses 
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The growing correlation length observed in supercooled liquids as their temperature is lowered has 
been studied with the aid of a single occupancy cell model. This model becomes more accurate as 
the density of the system is increased. One of its advantages is that it permits a simple mapping to a 
spin system and the effective spin Hamiltonian is easily obtained for smooth interparticle potentials. 
For a binary liquid mixture the effective spin Hamiltonian is in the universality class of the Ising spin 
glass in a field. No phase transition at finite temperatures is therefore expected and the correlation 
length will stay finite right down to zero temperature. For binary mixtures of hard disks and spheres 
we were not able to obtain the effective spin Hamiltonian analytically, but have done simulations 
to obtain its form. It again is in the universality class of the Ising spin glass in a field. However, in 
this case the effective field can be shown to go to zero at the density of maximum packing in the 
model, (which is close to that of random close packing), which means that the correlation length 
will diverge as the density approaches its maximum. The exponent v describing the divergence is 
related in d dimensions to the Ising spin glass domain wall energy exponent 6 via u — 2/(d — 2d). 

PACS numbers: 64.70Q-, 75.10.Nr, 64.70P- 
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I. INTRODUCTION 

One of the key concepts which has emerged in the last 
few years in the field of glasses is that of a growing corre- 
lation length scale £ [THS]. There are now many ways of 
defining and obtaining such a length scale: point-to-set 
0, patches [B], dynamics [U Uj etc. When it becomes 
large, they are probably all proportional to each other, 
as they are basically just a measure of the size of the co- 
operatively re-arranging regions in the liquid [5j . Simula- 
tions show that £ increases as the temperature decreases, 
or in the case of hard sphere and hard disk systems, as 
their density is increased. In this paper we report on our 
attempts to understand this growth, particularly in the 
context of hard disk systems in two dimensions but also 
for particles interacting with realistic potentials in any 
dimension. 

The leading theory for the growth of the correlation 
length has been that of the Random First-Order Transi- 
tion (RFOT) theory [THS]. In this theory the growth is 
driven by the decreasing configurational entropy or com- 
plexity PHI E] of the supercooled liquid as its temper- 
ature is decreased towards Tk, the Kauzmann tempera- 
ture [12] . In hard spheres there is a packing fraction 4>k 
at which the complexity apparently goes to zero, at least 
in the mean- field calculations of Refs. [TQ1HI]- At this 
density the correlation length diverges to infinity. How- 
ever, there are arguments that RFOT theory must be 
incorrect for systems in any finite dimension 13 . 

In this paper we shall try to understand the growth of 
the length scale £ not on the basis of RFOT theory but 
from lessons which have been learnt from studying in fi- 
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FIG. 1: (Color online) The hard disk system with single cell 
occupancy constraints. The square cells have grey outlines, 
and each cell can contain the center of only one disk (these 
are marked as black points). The outer edges of the disks do 
not interact with the cell walls, but only with the outer edges 
of other disks. Periodic boundary conditions have been used 
here and throughout this paper. 



nite dimensions the same p-spin models which inspired 
the RFOT theory. In Refs. [T3HT5] it has been shown 
that these models behave at low temperatures rather like 
an Ising spin glass in a field |TB] • Furthermore the corre- 
lation length grows as the temperature is decreased but 
saturates to a finite value at T = 0. It has also been 
argued that real glasses as well as p-spin models behave 
like Ising spin glasses in a field |T7]. This approach in- 
volved extensive use of the replica trick and is quite non- 
intuitive. It is one of the purposes of this paper to explain 
why, say, a binary mixture of hard spheres at high densi- 



ties will have features in common with Ising spin glasses 
in the presence of a field, but without the aid of the heavy 
machinery of replicas. 

To this end, we introduce in Sec. [il] the Single Occu- 
pancy Cell (SOC) model [18] [19] . In two dimensions it is 
a model in which the centers of the hard disks are each 
constrained to stay forever within a plaquette of a square 
lattice grid as in Fig. [I] (The generalization of this to 
higher dimensions is simple: in d = 3 one would use 
spheres whose centers are confined within the primitive 
cell of the simple cubic lattice). As the area of the disks 
is increased, the partition function of this constrained 
model becomes ever closer to that of the unconstrained 
model. This model with disks of the same size is not a 
glass: in fact it undergoes an Ising-like phase transition 
[20] to a state which is one of the two differently ori- 
entated slightly disordered crystal lattices shown in Fig. 
[3] In order to investigate glassy behavior we introduce in 
Sec. Ill a variant of the SOC model. This has two species 



of particles, A and B, present in equal numbers but ran- 
domly distributed over the plaquettes as indicated in Fig. 

m 

Fig. [2] also shows that the SOC model can be regarded 
as a spin model. An effective spin Hamiltonian is de- 
rived in Sec. |IV| for particles A and B which interact 
with a smooth potential V(r), e.g. the Lennard- Jones 
potential. For such potentials it is possible to calculate 
analytically a good leading order effective spin Hamilto- 
nian. The Hamiltonian is very familiar in the field of 
random magnetic systems; its vector spins have d com- 
ponents and interact with a d-component vector random 
field. The spin interactions are a mixture of exchange 
and pseudo-dipolar couplings and there are also single 
ion anisotropy terms. Because it is so well understood 
we shall just briefly outline in Sec. [V] the phases which 
can exist for the effective spin Hamiltonian. There are 
choices for the interatomic potentials for which the spin 
Hamiltonian is in the universality class of the Ising spin 
glass in a field and it this choice which is appropriate if 
one is interested in the properties of supercooled liquids 
or glasses [T7j . 

The Ising spin glass in a field does not have a phase 
transition in dimensions d < 6 [551 H3]> but the cor- 
relation length can become large as the temperature is 
reduced, provided the ratio h/J of the standard devia- 
tion h of the random field to the standard deviation J 
of the spin-spin coupling is small. In fact, we believe 
that for hard disks and spheres within the SOC model 
this ratio becomes zero as the packing fraction (density) 
approaches its maximum possible value </> max - The SOC 
model should therefore show features usually associated 
with ideal glass behavior in this limit. Much of the paper 
is devoted to investigating this intriguing possibility. 

Our analytical approach to the derivation of an effec- 
tive spin Hamiltonian does not extend to non-continuous 
potentials such as that appropriate to hard disks or 
spheres. In order to study them we have had to resort to 
simulations of the SOC model, in particular, event driven 




FIG. 2: (Color online) The spin mapping in action. Vectors 
are drawn from the center of each cell to the center of the 
disk occupying that cell - these vectors are the spins. The 
disks are then forgotten, and the system is treated as a spin 
system. Note that while there are two types of disk (large 
and small) there is only one type of spin after the mapping. 
The spin system contains quenched disorder, a consequence 
of the fact that each disk is constrained to remain in the cell 
to which it is first assigned. 



molecular dynamics. The details of this are described in 
Sec. EH 

We used the Lubachevsky-Stillinger (LS) algorithm 
[21] in Sec. VII to find some of the jammed states of the 
SOC model for hard disks. Its jammed states are simi- 
lar to those of the unconstrained model. At the densest 
packing possible, <^ mra , the state is jammed. We obtain 
an estimate of 4> ma , x from the largest value of the packing 
fraction 4>j of the jammed states which we have found 
in small systems, for which there is a chance that the LS 
algorithm might actually find the densest state. It is ac- 
tually very hard to do good simulations in the region of 



most interest, that is when 



because the con- 



straints introduced by the cell walls makes the dynamics 
even slower than that of the unconstrained system. In 
two dimensions, max turns out to be very close to es- 
timates of the glass close packing density (J)q, which is 
sometimes identified with the random close packing den- 
sity </> rcp [IT]. 

We study in Sec. 



VIII and Sec. IX correlation func- 



tions of the hard disk system in order to determine the 
variance of the random field h 2 and the variance of the 
spin-spin couplings J 2 . The physical reason for the pres- 
ence of a random field is also elucidated in Sec. IVIIII The 
form of the effective Hamiltonian is very similar to that 
obtained for smooth potentials in Sec. |IV| that is, it is a 
mixture of exchange and pseudo-dipolar couplings. Un- 
fortunately because of the difficulties associated with the 
long relaxation times as <j) — > </> max we cannot get good 
numerical estimates of how h and J vary with packing 
fraction in that limit. Fortunately we can provide an ar- 
gument in Sec. XII that shows h/J ~ (1 - 



One can use the droplet theory of spin glasses [2"5H2"T] 
to determine the growth of the correlation length £ from 
the ratio of h/J. According to the droplet picture, the 
correlation length £ can be estimated by equating the 
energy that can be gained from flipping the spins in a 
region of size £ in the random field, /i£ d//2 , to the domain 
wall energy cost of doing this, J£ , so 



£ 



which reduces for h/J^ (1 



* 



to 



(1 - (/)/(t) max ) 1 



d-26 



(1) 



(2) 



9 is the domain-wall exponent for Ising spin glasses in 
zero field. For d = 2, 9 w -0.287 J2S] so v « 0.78 while 
for d = 3, 9 « 0.24 (2§] and v w 0.79. Behavior of a 
power law kind is also expected in RFOT at a packing 
fraction <f>x < <f>Q. The value of v in that approach is 
dependent on whether or not "wetting" effects are con- 
sidered necessary [9]), but the wetting form v = 2/d is 
not very different from that of Eq. (p ) in two and three 
dimensions due to the fact that in these dimensions 9 
is small. However, in our approach, we have not seen 



any evidence for the ideal glass transition at 4>k- For 
us the divergence of the correlation length is associated 
with glass close packing and jamming. 

Finally in Sec. |XIII| we discuss the key question; which 
features of supercooled liquids and glasses can the SOC 
model be expected to describe correctly? It is argued 
that the SOC model should be good for understanding 
some of the phenomena which exist on time scales less 
than the alpha relaxation time, as the caging of the par- 
ticles on time scales less than the alpha relaxation time 
is mimicked by the trapping of the particles in the cells 
in the SOC model. The dynamical correlation length is 
extracted from the properties of correlations at the al- 
pha timescale so we expect that the SOC model should 
at least give v correctly. 



II. THE SINGLE OCCUPANCY CELL MODEL 

Cell occupancy models have a long history in the study 
of phase transitions in fluids and liquids [HI [19] ■ In the 
past, they have been used to calculate the equation of 
state of hard spheres at high density [30] or to place 
bounds on derivatives of the free energy |31| or entropy 
[I9l . The system is divided into cells of a chosen geom- 
etry and a constraint is applied which fixes the number 
of particle centers found in each cell. We focus on the 
single occupancy cell (SOC) model, where each cell can 
contain at most one particle. We work in two dimensions, 
although the model is easily generalized to higher dimen- 
sions. Fig. [I] shows a hard disk fluid with a single cell 
occupancy constraint using square cells. The constraints 
mean that disk centers interact only with cell walls and 
disk surfaces interact only with other disk surfaces. 

We note that it might be possible to realize the SOC 
system experimentally, at least in two dimensions. The 
square cells could be produced by a wire grid, and a post 
could be attached at the center of each disk so that while 
the circumference of the disk can pass under the wire 
grid, the post at the center cannot. 

SOC models are useful to us because they make the 
introduction of a spin representation of the problem 
straightforward. A disadvantage of using the cell con- 
straint is that at low packing fractions the behaviour of 
the system deviates significantly from the behaviour of 
the unconstrained system. At low packing fractions most 
of the collisions will be between disks and cell walls, so 
the cell geometry dominates. As the packing fraction is 
increased, more collisions occur between disks and close 
to jamming, almost all collisions will be between disks. 
The closer the packing fraction is to <^ ma x, the better an 
approximation the constrained model becomes to the un- 
constrained model as the cell walls no longer dominate 
the dynamics. 

Another pecularity of SOC models is the appearance of 
singularities in thermodynamic properties. This occurs 
because of how the constraints limit the size of clusters 
that can form. Without constraints, it is possible to find 



all particles forming a single cluster at all packing frac- 
tions. This is not possible in the constrained system. As 
shown by Hoover and Alder for a one-dimensional hard 
rod SOC modcl[18 , it is only possible to form clusters 
of a certain size above a certain packing fraction. For 
example at very low packing fractions, the constraints 
mean that clusters can only contain at most two particles. 
As the packing fraction is increased clusters can contain 
three then four particles. At the packing fractions where 
it becomes possible for larger clusters to form, the par- 
tition function changes its analytic form and this means 
that discontinuities appear in thermodynamic quantities 
such as d 2 P/dV 2 . These packing fractions get closer 
together approaching <f> max , and the discontinuities de- 
crease in size, meaning that the shortcomings become 
less important. We expect similar behavior in two di- 
mensions, but as each disk has more nearest neighbors 
the effect will be smaller. In any case, it has not been 
noticeable in our simulation results. 

The model with disks of the same size, as in Fig. [T] 
is not a glass. Without constraints the largest possi- 
ble value of the packing fraction <p occurs when a tri- 
angular lattice with all disks touching there neighbours 
is formed; so </> max = tt/2V3 « 0.9069. In the SOC 
version of the model, the constraints mean an exact tri- 
angular lattice cannot form, so we hnd (f> max ~ 0.88. 
When </> = <fi c w 0.77 there is a phase transition to a 
slightly disordered crystal which is orientated in one of 
two possible directions as in Fig. [3] The critical ex- 
ponents of this transition are expected to be those of 
the two-dimensional Ising model because of this two-fold 
degeneracy of the orientation of the slightly disordered 
crystal lattices [20] , 




III. SINGLE OCCUPANCY CELL MODELS FOR 
MODELLING GLASSES 



FIG. 3: (Color online) A snapshot of two 16 x 16 hard disk 
systems in the SOC model at (f> = 0.83, each displaying one 
of the two possible defected crystal orientations. 



To make a glassy model, we introduce two different 
sizes of disk. The binary disk system consists of hard 
disks of two species (A and B), where the SOC con- 
straints have been applied, as in Fig. [2j The species 
of disk have different radii ex a and <tb , where the size ra- 
tio i?AB = cta/ctb is held fixed as the packing fraction is 
altered. The packing fraction is given by 



cf> = n(a 2 A + a|)/2, 



(3) 



where the side of the square plaquette has been taken to 
be of unit length. We set Rab = 1.0/1.4 - this is a well 
explored choice [33H35] . There are equal numbers of each 
species (Na = Afe = N/2), and each cell contains a disk 
of species A or B with equal probability. 

When Rab -^ 1, the system undergoes crystallization 
to one of the two disordered crystal states shown in Fig. 
|3]but with substitutional disorder. Disks of species A and 
B will be distributed at random throughout the defected 
crystal. 



For Rab = 1.0/1.4 without constraints, the densest 
state is a phase separated crystal where the two species 
form separate triangular crystals. Although this state is 
very stable, it takes such a long time to form that it is 
rarely reproduced in simulations. This makes the system 
a good model glass former. Recent work has shown that 
phase separation may be achieved on simulational time 
scales in some three dimensional binary systems [37] • It 
may be that some nacscent phase separation could be 
driving behaviour normally identified as glassy (slow dy- 
namics, dynamic heterogeneity and growing correlation 
lengths). For an example, see Ref. [55] , 

With the introduction of the single occupany con- 
straints, phase separation can no longer occur as fixing 
the species of the disk in each cell fixes the local com- 
position of the hard disk fluid. We chose to distribute 
the species across the cells with equal probabilities. This 
mimics what would happen if a well mixed fluid at low 
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FIG. 4: (Color online) Mapping of a disk to a spin - the 
spin Si is defined by the equation Xi = Ri + Si where Xi is the 
position vector of disk i and Ri is the position vector of the 
center of cell i. 



packing fraction was rapidly 'quenched' to a higher pack- 
ing fraction without allowing the disks to phase separate. 
As phase separation is prevented by the cell constraints, 
this means that glassy behavior can be investigated in a 
fully equilibrated model; there are no concerns that if one 
runs the simulation for longer there will eventually be 
phase separation. 



IV. THE SPIN HAMILTONIAN 

Our main reason for studying the SOC model is that 
it makes mapping to a spin system easy. This is acheived 
by drawing a vector from the center of each cell to the 
center of the disk that occupies that cell. The cells are 
labelled i, where i — 1,2, ... ,N and the spin Si is defined 
as 



Xi -**i i $i: 



(4) 



where Xi is the position vector of the disk i and Ri is the 
position vector for the center of cell i. This mapping is 
illustrated in Figs. [2] and [4} 

The system may now be is analyzed as if it were a spin 
system. Unlike the particles, the spins are equivalent and 
all details of the interaction between them are found in 
the terms and couplings of the spin Hamiltonian. If we 
write Si = (xi,yi), which is appropriate for d = 2, then 
in the usual XY model \si\ = 1, but in the SOC spin 
model, Xi and yi take values which keep the disk in the 
ith plaquette. 

Let us suppose that the particles in the cells interact 
with each other through the potential V(r), where r is 
the interparticle separation. The hard disk problem is a 
special case of this potential where V(r) = oo if r is less 
than the sum of the radii of the two disks and is otherwise 
zero. If we have a binary mixture of two types of particles 



A and B, V(r) will be a shorthand for Vij(r), where i and 
j encode the species of the interacting particles. We shall 
now proceed to derive the effective Hamiltonian in terms 
of the spin variables Si . 

Using the notation in Fig. [4j the distance r between a 
particle in cell i and one in cell j is 

r = \Si - x 3 \ = \Ri + s t - Rj - Sj | . (5) 

To second order in the spin variables, 

V(r) = V(R) + V'(R)R-(s i ^s :J )/R 

+ [V"(R) - V'(R)/R] \R ■ (^ - Sj)] 2 /(2R 2 ) 



+ V'(R) [^ - s 3 f /(2R) + 



(6) 



where R = liL- — R« 



Ri 



The partition function Z of the SOC model is 



f N 

/ JJ dxidyi exp[— 0H] 



(7) 



where the integration over Xi,y% covers the area of the 
ith plaquette. The Hamiltonian % is, to second order in 
the spin displacements, of the form (up to constants) 



u = -J2K^~ 



1 



E 



(8) 



where the sums over /j, and v run from 1 to d and in 
d = 2, ^ = (xi,yi). The fields hf are given by 



K 






)K/Ri. 



(9) 



If all the particles are identical, the "field" term h^ is 
identically zero. However, if we have a binary mixture 
of two types of particles A and B such that Vaa, Vbb, 
and Vab all differ, then the field term h^ is non-zero and 
time-reversal invariance is broken. 

We can calculate the average of h^ when the average 
is taken over the various possibilities allowed by the se- 
lected disk distribution. We will consider just nearest- 
neighbor interactions to illustrate how the calculations 
can proceed, and the case v — x. Only the sites to the 
right and left of the site i contribute to the sum in Eq. 
(p)]). At each of these sites there can be an A or a B disk 
(with equal probability) and at the site i there is an equal 
probability of the disk being A or B. Summing over the 
various possibilities one finds h% = 0. The distribution of 
the random field components is such that h^h" = h 2 8^ w , 
where 

h 2 = \[{V AA - V' AB f + (V BB V BB n (10) 

The various derivative are calculated at the nearest- 
neighbor distance. Note that if V A a — Vbb = Vab, 
then h = 0, as expected. 



The quadratic term in Eq. Q takes the form for i ^ j 



H eff (s i ,s J ) = ^-> 



Si ■ s 3 - [1 - 



Ri 

RjjV"{Rii) yfi -U6 -N 

v(%) ](^- s «)(%-^-) 



> (11) 



where the unit vector 12^ is Rij/Rij. Note if the inter- 
action V(r) = — A/r n , Eq. ( |Tl"| ) reduces to 

« J r ^ ^* - 

H eff (s t , Sj) = -—^+2 ^ ■ g 3 - ( U + 2 )( R *3 ■ Si){R t j ■ Sj) 
&ij L 

which for n = 1 is the familiar dipole-dipole cou- 
pling interaction. For other non-power law potentials 
H e ff(si, Sj) can be regarded as a mixture of the exchange 
interaction with pseudo-dipolar couplings. 
When i y^ j, D^ is of the form 
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rjM" AX _ R P^ U v 

^ij ~~ sllj< J fll> JD l jJX i j£l i j 



A, 



V'(Ry) 
7? * 



(12) 



(13) 



FIG. 5: (Color online) Cells 1 to 20 contain disks which can 
possible collide with the blue disk when the hard disk fluid is 
at high packing fraction. 



and 



B, 



V(Rij)r« _ RtjV"(Rij)-\ 



Ri 



V'(Rij) 



(14) 



For i = j there are single-ion anisotopy terms with coef- 
ficients 



D 



fits 



= -J2 D 



ij ' 



(15) 



J¥* 



For smooth potentials like the Lennard- Jones potential, 
the h^ and D?? can therefore be directly calculated. The 
configurational average and variance of D?f due to the 
quenched random distribution of A or B particles in the 
plaquettes can be obtained by the method used to obtain 
Eq. (10); the expressions are complicated. 



The only approximation which arises from the use of 
the Hamiltonian in Eq. ([8]) is the truncation to second 
order in sf. The hope is that this truncation does not 
alter the "universality class" associated with the phase 
transitions of the spin system. Of course, further terms 
could be included if required. 

When Eq. (15) is used to fix the single-site terms, 



the Hamiltonian will still have in its quadratic terms the 
translational invariance of Eq. ffib . Similarly if 



W 



/ sCjjRjj, 



(16) 



that will ensure translational invariance in the linear term 
in Eq. (p|. The quantities A^j, B t j and CV,- thus specify 
an effective spin Hamiltonian for our problem. 



In Fig. [5] the positions are shown of the plaquettes 
whose associated disk can interact with the disk in the 
central plaquette when the packing fraction is high. The 
number of such disks is surprisingly large; 20. At smaller 
packing fractions the number is reduced to 8. (In three 
dimensions the number at large packing fractions is 80). 
Now for the blue disk to interact with the disk in plaque- 
tte 1, the disks in 2, 5, 6, and 10 must be occupying only 
a restricted portion of their plaquettes. A complicated 
many-spin set of terms in the effective spin Hamiltonian 
is needed to describe this feature. It is clear that keep- 
ing for example only nearest-neighbor spin-spin interac- 
tions does not contain the physics of the increase in the 
effective number of interacting spins as the packing frac- 
tion increases. Truncating the effective Hamiltonian to 
just binary spin interactions may also fail to capture the 
properties successfully modelled by p-spin models such 
as the dynamic transition. In this paper, our main con- 
cern is the behavior of glasses at temperatures below the 
dynamic transition temperature or at densities above <f>d, 
the packing fraction associated with the (avoided) dy- 
namic transition (see Sec. XIII ) and binary spin interac- 
tions are quite sufficient to capture the Ising spin glass 
behavior which prevails there. An investigation as to 
whether the considerable range of the spin interactions 
can explain the utility of mean-field ideas in glasses [TT] 
is being carried out 38 . 

For hard disks and spheres the potential V(r) is infi- 
nite when r is such that they overlap, and zero otherwise. 
Such a potential makes V'(r) zero except at the contact 
distance where it is infinite. As a consequence the ex- 
pansion used in Eq. (|6| is not useful. For hard disks 



and spheres we shall still use quantities like h? and D^ , 
but instead of deriving them from the potential we will 
obtain their values as fitting parameters chosen to repro- 
duce measured correlations (like (sf) = (x/) etc.) rather 
in the spirit of Ref. [35]. This is done in Sec. IX 

We have already noted that when Rab — > 1, the sys- 
tem will undergo crystallization to one of the two disor- 
dered crystal states in Fig. [3] but with substitutional dis- 
order. The disks of species A and B will be distributed 
at random throughout the defected crystal. However, 
when Rab is close to 1 there will be effectively random 
fields arising from the small differences in the A and B 
particles. We suspect that this changes the transition to 
the disordered crystalline state to that of the random- 
field Ising universality class. We shall suppose from now 
on that Rab is sufficiently different from unity that this 
crystal-like transition no longer arises and that only glass 
ordering behavior (i.e. spin-glass ordering in the spin 
mapping) need be considered. 



V. SPIN GLASS BEHAVIOR 

In this section we shall discuss the properties of a spin 
Hamiltonian like that in Eq. ([8]). For any smooth poten- 
tial, the hf and Df/' l/ can be directly calculated from the 
potential. These expressions will be renormalized by the 
effects of multi-spin interactions neglected in Eq. pb, 
but hopefully they provide a good first approximation. 
For hard spheres or disks they are parameters obtained 
by fitting to the measured correlation functions (see Sec. 

ix}- 

For the binary SOC model, the Df" between sites i 
and j will depend on whether the particles in the plaque- 
ttes are A or B particles. As the particles can never es- 
cape from their cells, there is quenched disorder present. 
One can obtain the probability distribution function of 
the W and Dff? and obtain their mean and variance. 
Rather than do this, (which is rather cumbersome and 
uninformative), we will just outline some of the possibili- 
ties which might arise. What actually happens for a given 
set of potentials Vaa, Vab, Vbb requires explicit calcula- 
tions and simulations and the number of possible phases 
is large. To limit the discussion it is useful to recall the 
underlying system: disks (or spheres) whose centers are 
trapped in the squares (cubes) of a square (simple cu- 
bic) lattice. Ferromagnetic ordering in the spin system 
would correspond to a crystallization of the disk centers 
into a square lattice of the same periodicity as that of 
the plaquettes. This will not happen if one uses an ap- 
propriate binary mixture for modelling glasses and so we 
will discount the possibility of a transition to a ferromag- 
netic state and just concentrate on situations which are 
spin-glass like, i.e. those where the standard deviation of 
the couplings D^ dominates their mean values. We will 
therefore not be discussing the type of ordering shown in 
Fig. k^ (for the monatomic system) where there is clear 
crystal order present: Glass behavior is not associated 



with any kind of long-range crystalline order. 

The spins in the system are (i-component spins so that 
one might have thought that any spin glass phase in 
this system would be in the universality class of the d- 
componcnt vector spin glass. However, it was shown 
a long time ago [?D] that in the presence of pseudo- 
dipolar-like terms, the transition to the spin glass phase 
is changed from one in the d- vector spin glass universality 
class to one in the Ising spin glass universality class. 

The random field terms hf have a dramatic effect on 
the nature of the spin glass. At mean- field level and for 
dimensions d > 6 a d-component random field present in 
a d-component vector spin glass produces a phase transi- 
tion - the de Almeida-Thouless transition [2T] - which is 
in the same universality class as an Ising model in a field 
[2"2l HI] . The presence of the pseudo-dipolar terms just 
reinforces the Ising nature of this transition. For d < 6 
the spin glass transition is removed by the presence of 
the random field [22| [23]. 

The spin-glass correlation length, which is equivalent 
to the point-to-set length scale, can still become large for 
d < 6 if the ratio h/J is small. (J is a measure of the 
standard deviation associated with the D^). According 



to the droplet picture [25-27 the correlation length £ de- 
pends on this ratio as in Eq. M. This is the correlation 
length appropriate to T — 0. As a function of tempera- 
ture the correlation length is small at high temperatures 
and grows to this value in the limit when T — > 0. We ex- 
pect that £ might become large for real fragile glasses at 
low temperatures. However, on this picture £ will never 
become infinite unless h/J goes to zero. We suspect that 
this never happens for smooth potentials. In other words, 
for such potentials no diverging length scale is expected 
in d < 6. 

Note that if we had used the mean-field approximation 
to determine the properties of the spin system, we would 
have found a phase transition, the de Almeida-Thouless 
transition, at a finite temperature provided the ratio h/J 
is not too large. We would have then been tempted to 
identify this transition with the ideal glass transition. 
However, it is our belief that the AT transition does not 
occur for dimensions d < 6 [22J [23] . 

One might further wonder whether the multi-spin "p- 
spin" interactions which were alluded to in the discussion 
of Fig. [5] might make a transition to a one-step replica 
symmetry broken state possible. We have neglected them 
in our discussion. This is the scenario envisaged in the 
RFOT and is the origin of the ideal glass transition. We 
do not think such a transition can exist outside the mean- 
field approximation, that is, in finite dimensions, where 
the one step replica symmetry broken state is unstable 
against the thermal excitation of large droplets [32] ■ 

For binary mixtures of hard spheres and disks, a mech- 
anism might exist to drive the ratio h/J to zero. In the 
SOC model there is a maximum packing fraction for hard 
disks or spheres. For our binary mixture of hard disks, 
this value is estimated in Sec. |VII| Its value (/> max is very 
similar to 4> rcp of the unconstrained model and in both 



models at these densities, the pressure is infinite. We 
shall present numerical evidence and arguments in Sec. 
IX that the ratio h/J ~ (1 — </>/</> max ), so that in our ver- 
sion of the SOC model, the correlation length £ diverges 
to infinity according to Eq. (TH). In other words there 
are features of a glass transition in the hard disk SOC 
model as <p ~ * </>max! i n that there is a diverging correla- 
tion length in spin-glass-like correlation functions. The 
rest of this paper is devoted to the study of this behavior. 
To acquire data to determine </> max and to obtain esti- 
mates of h and J, it is necessary to perform simulations 
of the SOC hard disk system. In the next section, our 
simulation methods for hard disks are outlined. 



VI. EVENT DRIVEN MOLECULAR 

DYNAMICS AND THE 

LUBACHEVSKY-STILLINGER ALGORITHM 

To simulate the hard disk system, we use event driven 
molecular dynamics following the method described by 
Lubachevsky [44] . This is an efficient way to perform 
simulations of hard disk systems. We will not describe 
the method in full here, but the basic principle involved 
is to keep a list of the next collision each particle will 
be involved in ordered by time. Time is moved forward 
by jumping to the collision that occurs next, and then 
recalculating the list in light of the new velocities and 
positions the colliding particles now have. The speed 
of the simulation is further boosted by the fact the cell 
constraints restrict the particles that can possibly collide. 

We generate configurations at a particular packing 
fraction by first placing particles randomly in each cell 
in such a way that each cell is equally likely to contain 
a particle of either species. The particles start with zero 
radius (so there is no possibility of overlap) and at time 
t have radius <7j(t) = Tit where the i denotes the species 
of the particle in question. The growth rate I\ is set to 
be small to allow the disks to remain in equilibrium as 
the simulation progresses. We use T ~ 10~ 4 . Each disk 
is given a random velocity so that \vi\ = 1 

The disks are allowed to collide and grow until the 
system reaches the desired packing fraction. Then the 
disk radii are set to be constant and measurements may 
be made. 

To generate jammed states, we make use of the 
Lubachesky-Stillinger algorithm [53] . We begin the sim- 
ulation as described above, but in this case the growth 
rate of the disks is not set to zero at any time. As the 
simulation proceeds, collisions (events) become separated 
by smaller and smaller time intervals and the simula- 
tion will become slower. If Stuu is the time between 
disk-disk collisions, l/Stno — > oo as the simulation pro- 
ceeds. This is equivalent to a divergence in the pressure. 
The Lubachevsky-Stillinger algorithm works by choosing 
a value of St* DD below which collisions are close enough 
together that the system has effectively jammed. Here 
we use St 



ulation yields a range of jammed configurations, with a 
distribution of jammed packing fractions </>j. 



VII. DETERMINING THE MAXIMUM 

PACKING FRACTION IN THE SOC BINARY 

MODEL 

In this section we shall estimate the largest packing 
fraction (f> max for our binary hard disk system. It is as 
~> '/'max that we expect the correlation length to di- 
verge, so max is like the critical temperature of the sys- 
tem. 

As already discussed, at low packing fractions the be- 
havior of the constrained fluid is very different from that 
of the unconstrained fluid, becoming closer to it as the 
packing fraction is increased. At some packing fraction 
the system will jam. In a jammed state for the uncon- 
strained system, the disks are held in place by their z 
neighbors, (except for a few rattlers), where z — 2d - 
the so-called isotatic condition [35]. In the SOC model a 
disk can be jammed when its center is pinned against a 
plaquette wall. 

Jammed states were obtained for a range of system 
sizes using the Lubachevsky-Stillinger algorithm [23] de- 
scribed in Sec. VI A plot of the values of <j>j values for 



DD 



10 with {\vi\) — 1. Repeating the sim- 



the binary hard disk system can be seen in Fig. [6] The 
plot does not show the complete range of jammed states 
possible in the L x L, but a subset obtained from several 
runs of the algorithm. 

There is a fall from </> max ~ 0.855, when L = 5, 
then (p max (L) settles around <f) = 0.835 before falling off 
slightly when L > 30. This fall off is due to inefficien- 
cies in the simulation for generating jammed states of 
the highest packing densities. The value <j> max sw 0.835 
is quite close to the value of the packing fraction at 
which unconstrained binary disk systems of this type 
jam, (f>j sa 0.84, using the protocol studied in Ref. [33] . 
In other words, it is close to the numbers quoted for "ran- 
dom close packing" in two dimensional systems. 

The most significant change from the unconstrained 
binary fluid is the presence of a well-defined maximum 
jamming density max . In the unconstrained fluid when 
a jammed configuration has been acquired, one can al- 
ways imagine creating a denser state by rearranging a few 
of the particles into a region with more local crystalline 
order. This will create a small amount of free volume 
which will allow further arrangements to be made. If 
this programme is continued, the final point is a com- 
pletely crystalline configuration. A continuum of states 
at packing fractions between <pj and 4> C rystai can be con- 
structed by this method (although there is no guarantee 
that they will be stable). This makes defining a dens- 
est non- crystalline state problematic. However, since the 
cell constraints do not allow the composition of the fluid 
to be altered, this programme cannot be followed in the 
SOC model and there is indeed a well-defined maximum 
density. This maximum density will depend on the par- 
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FIG. 6: The jammed packing fractions for the SOC binary 
hard disk fluid found in simulations from different random 
starts at various linear dimensions L. 



Imagine a binary fluid at high packing fraction, fo- 
cussing on one single disk. At any given time there will 
be a variety of moves the disk will be able to make. Most 
will be short and rapid (the typical behavior of a caged 
particle), but some may be part of large rearrangements 
that will allow the structure of the fluid to relax and 
change. It is reasonable to assume that the cell con- 
straints will block a lot of these movements (simply be- 
cause the walls of the cell will intercept the paths the 
disk wants to take) , and they are more likely to interfere 
with the longer paths. Thus with the cell constraints 
in place, it is expected that the dynamics of the system 
will become slower. Lots of local rattling will be allowed, 
but the system will have to wait for longer before large, 
co-operative movements that allow structural rearrange- 
ments take place. 



VIII. ORIGIN OF THE RANDOM FIELD FOR 
HARD DISKS 



ticular realization of the distribution of large and small 
disks over the cells, but it is probably a self-averaging 
quantity. 

In the unconstrained model each protocol for produc- 
ing jammed states produces states with a characteristic 
value of 4>j as N — > oo. The Lubachevsky-Stillinger al- 
gorithm used in this paper produces, in the SOC model, 
states of a characteristic <j)j, which will not in general 
include the states at ^ max , except possibly at small val- 
ues of N. Within the SOC, different protocols will also 
produce different values for <j)j. Protocols which produce 
jammed states whose <f>j is close to max are produc- 
ing jammed states closer to those in the unconstrained 
model. As a consequence, we are expecting that for the 
densest jammed states, nearly all the disks will be touch- 
ing 4 other disks in the jammed state and very few, if 
any, will be jammed because their center is touching a 
plaquette wall. In principle, but probably not in prac- 
tice, one could obtain estimates of max by calculating 
the pressure P in a fully equilibrated system and deter- 
mining (f> max by fitting to 



PV 



Nk B T </>max(l - </>/^max) 



(17) 



which becomes exact as <j) — > </> max ,30 . The problem 
with using this procedure is that it is very hard to equi- 
librate the system at packing fractions close to </> max - 

The cell constraints affect the dynamics of the system. 
This is a key concern as it affects how quickly the system 
can be equilibrated and hence the quality of simulations 
which can be done. The system is clearly glassy - simula- 
tions performed on systems with packing fractions above 
4> = 0.75 become noticeably slow, while approaching the 
maximum packing fraction of around </> s=s 0.835 makes 
good measurements extremely hard. The presence of the 
cell constraints makes the dynamics even slower than that 
of the unconstrained system. 



We have already remarked that when all the particles 
are identical the field term /if in Eq. (|8| is zero. For 
binary mixtures it is non-zero and this makes the expec- 
tion of the local magnetization (sf ) also non-zero. This 
is easily understood from Fig. [7] 

Zero local magnetisation means that a disk spends its 
time symmetrically distributed over its cell. With this 
in mind it is easy to see why the local magnetization 
is finite at all (f> in the binary system. When the pack- 
ing fraction is very low the disks rattle backwards and 
forwards in their cells, rarely colliding with each other. 
The finite local magnetization is caused by having disks 
of different sizes on either side of the central disk. Say 
there is a large disk to the right, and a small disk to the 
left (as shown in Fig. uh. The neighboring disks will in- 
trude into the cell. When their sizes are different they 
can intrude by different amounts. In the case just de- 
scribed the central disk will spend more time on the left 
hand side of the cell as there is more free volume there. 
As the packing fraction is increased, there is more intru- 
sion by the neighboring disks and the deviation from the 
center of the cell becomes larger This means that the lo- 
cal magnetization gets larger. This suggests that there 
should be three different types of behavior for the local 
magnetization: large disk to the left and small disk to 
the right ((sf) > 0), large disk to the right and small 
disk to the left ((sf ) < 0) and lastly disks of the same 
size on each side ((sf) ~ 0). In Fig. M the components 
of the local magnetisation split into these three groups. 
At higher packing fractions, the groups blur into one due 
to interactions between increasing numbers of disks, and 
the components are randomly distributed about zero. 

In the spin interpretation, a finite local magnetization 
randomly distributed about zero implies the presence of 
a local random field hi interacting with each spin through 
a term of the form — J^. hi.si. The expectation value of 
the total magnetization (M) = 0, where M = J^i sl/N , 
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FIG. 7: (Color online) Origin of non-zero local field and 
hence magnetization (sf) in the binary disk system. Disks of 
different size on either side of the central disk bias it in one 
direction - in this case the larger disk on the right forces the 
central disk to spend more time on the left of its cell. 




FIG. 8: (Color online) All components of local magnetization 
(sf } for the spin system derived from the binary hard disk 
fluid under the SOC at low packing fraction for a 4 x 4 system. 
Note the three distinct bands into which they fall. 
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FIG. 9: (Color online) Averaged spin-spin correlations (sf sj) 
and (sfs?) between spin i and its nearest neighbors and next- 
nearest neighbors, averaged over sites i and disk realizations. 
Note that sets of points with the same color all lie on top of 
each other. 




for all packing fractions. This suggests that the /if must 
be evenly distributed around zero. The source of the 
random nature of the field is the random distribution of 
the species of disk over the cells, since this affects the 
local magnetization at all packing fractions through the 
mechanism described above. This field will be discussed 
again in the following sections. 



IX. SPIN-SPIN CORRELATIONS 



We now study correlations of the form (s^s^) where 
i and j index the lattice sites and v and /i label the 
x and y components of the spins. For a spin in cell i, 



we can calculate 



Isfsf) 



ksU)) 



ktfs)) 



and {sis*) for 



nearest neighbors (the spins north, south, east and west 
of spin i) and next-nearest neighbors (the spins north- 
east, south-east, south-west and north-west of the spin 
i). We are interested in using these correlations as a 
guide to the effective interaction between the hard disks. 



Our studies suggest that the effective spin interactions 
follow closely the form expected in Sec. |1V| the effective 
spin Hamiltonian is well-approximated by Eqs. ([8]) and 



There are many different local environments a disk can 
experience. We have therefore studied the average of 
these correlations, defined as follows. We have calculated 
for each site i its spin's correlation with its neighbors at 
i+S, where S is a label running over the N,W,E,S nearest 
neighbors and NW, SW, SE, and NE next-nearest neigh- 
bors (i.e. we calculate (sfs" +s )). The site averages of 
these correlation functions were also calculated and the 
results are shown in Figs. [9] and [10] for N = 256. In 
principle there is no need to do an average over disk real- 
izations as the site averages are self- averaging quantities. 

There are some notable features visible in these Fig- 
ures. The correlations are seen to grow as the packing 
fraction increases, suggesting that the coupling between 
spins increases in strength with packing fraction. Study- 
ing Fig. |9j the strongest correlations are seen to be those 
with the East and West spins for (sfs 3 -) and with the 
North and South spins for {s^s^}. It makes sense that (for 
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FIG. 11: Overlap q measured for a range of packing fractions 
in a system with N = 256. 



The overlap measured for a spin glass in a field is finite 
at all values of the temperature, growing larger as T — > 
0. This happens because the (random) fields bias the 
orientations of the spins. 



X. EFFECTIVE SPIN HAMILTONIAN FOR 
HARD DISKS 



FIG. 10: (Color online) Averaged spin-spin correlations 
(sfs^) and (sfsj} between spin i and its nearest neighbors 
and next-nearest neighbors. 



example) when a disk is moved to the East, its neighbors 
to the East and West should also move in that direction. 
This will generate the large correlations observed when 
(s^s'j) is measured with the spins to the East and West 



of the central spin. Studying Fig. 10 it is clear that for 



</\ 



{sfSj) and (s^Sj), the North, South, East and West cor- 
relations are all zero while the others are small but finite. 
This confirms the presence of pseudo-dipolar interaction 
terms in the effective Hamiltonian and is compatible with 
a Hamiltonian like that of Eqs. (l8|) and (12 1. 



We have also determined the Edwards- Anderson order 
(overlap) parameter, defined as 



<1 = 



N 



oV\*\ 



(18) 



where the square brackets [. . .] av mean an average over 
the quenched disorder in the system (here the species of 
particle that each cell contains) . The overlap is a measure 



of the amorphous or glass order in the system. In Fig. 11 
it can be seen that the overlap increases as the packing 
fraction is increased towards its maximum possible value. 
It is always non-zero even at small packing fractions. 



We will try to understand the correlations studied in 
Sec. HXI with the aid of an effective Hamiltonian like that 
in Eqs. pi) and ( [12] ), but for simplicity we ignore all 
couplings except those between nearest neighbors. This 
is a poor approximation at large packing fractions, but 
is better for low packing fractions. We will also work to 
lowest non-trivial order for each quantity studied. 

A weak coupling expansion can be made which allows 
fitting of Aij and By from the simulation results. Unfor- 
tunately, as this is a weak coupling approximation (i.e. it 
is valid when Aij and By are small) it cannot be used to 
accurately measure them in the region of most interest, 
<t> — > 0max, as there they become large. 

The correlation < sjV' > is calculated using: 
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v M e ~0H eff 
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n ^ ( ig ) 



where the spins components s^ are integrated over the 
fcth cell, which has unit side length. The partition func- 
tion is 



-1/2 /.1/2 



/ / e-P*'** TT rf ; 

J -1/2 J -1/2 



Sfc- 



(20) 



The integrals can be performed by first Taylor expanding 
the exponential, and then integrating to give the correla- 
tion in terms of Ay, and J5y and some simpler averages. 
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FIG. 12: (Color online) Fitting of the averaged 0Aij and 
j3Bij for a binary hard disk system under SOC with N = 256 
particles through spin-spin correlations using Eq. ( |21[ ). The 
four curves for each quantity come from looking at the North, 
South, East and West directions. Local environments have 
been averaged over leading to a strong symmetry between 
the directions. 
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and expanding 
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Again one can replace ((s^) 2 ) by 1/12. 
Thus the variance h 2 , denned as 



(23) 
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72q. (24) 



Eq. ( 24 ) shows that the variance of the random field will 



increase with packing fraction just like q does, at least 
when q is small, (see Fig. 11). The equation will not 



hold at high packing fractions, where we actually expect 
{/3h) 2 to diverge but q must always remain less than i. 
(This inequality arises because q cannot exceed the value 
it would have if all the disks were simultaneously at the 
corners of plaquettes). 



XI. THE CORRELATION LENGTHS 



On perforating the Taylor expansion we find 

«l/2 /-1/2 
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and performing the integration yields 
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(21) 



To the order we are working («) 2 ) = ((sf) 2 ) « 1/12. 

Because we have measured {s"s^) we can use these 
measurements to determine f3Aij and j3Bij for each bond 
(nearest-neigbor pair). The values of (3Aij (and fiBij) 
have a distribution, with a mean and a standard devi- 
ation. The standard deviation is important as it is the 
randomness of the effective couplings which is encoded 
in the standard deviation which can be the cause of spin 
glass behavior if it is sufficiently large compared to the 
means of the couplings. In Fig. [12] we have plotted the 
averages of j3Aij and /3Bij as a function of the packing 
fraction <fi. 

Using the same approximation for the effective Hamil- 
tonian we can determine the variance of the random field 
hi from our results for q. 
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(22) 



The quantity of most interest is the spin glass correla- 
tion length as it should be the glass correlation length. 
We shall determine it via the spin-glass susceptibility. 
First the cumulant xtj = ( s i s j) ~ ( s i)( s j) is obtained. 
This measures fluctuations in the correlations between 
the fj, and v components of the spins i and j. The spin- 
glass wave-vector dependent susceptibility is [U] 
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* ( fc ) = Tfl^L l lWZrUexp(ik.Ri j ), 
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where Rij is the vector connecting lattice sites i and j. 
From x SG (k) the spin glass correlation length £ SG can 
be calculated using the formula 46 
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(26) 



where k m i n is the minimum non-zero wave- vector k m i n = 
(2tt/L,0). 

Additionally, a ferromagnetic correlation length can be 
calculated and compared to the spin glass length to see 
which kind of correlations are dominating the system. 
A ferromagnetic wave-vector dependent susceptibility is 
defined: 
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This is similar to the spin glass susceptibility, but xtj is 
not squared. From this, a ferromagnetic length scale £ F 
can be calculated: 



2sin(|fc mi „|/2) 
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(28) 
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FIG. 13: (Color online) Spin glass correlation length (in red) 
and ferromagnetic correlation length (in green) for the spin 
system derived from the binary hard disk fluid. Calculated 
for a system with N — 256 disks. 



Fig. [13] shows these length-scales plotted together. It 
is interesting to study both cases, since with the map- 
ping to a spin system there is not yet an a priori way 
of predicting the properties of the spin system. There 
exists another mapping of the structural glass to a spin 
system by Stevenson et al. [47] ■ This method is similar 
to the mapping of Moore and Yeo [T7] , in that it requires 
a replication of the structural glass Hamiltonian but it re- 
sults in a random bond Ising model in random field. The 



for the system to equilibrate and the susceptibility to 
reach its correct level. 



XII. SCALING OF THE EFFECTIVE 
COUPLINGS NEAR d> mBX 



As 4> — > </> max it is just not possible to equilibrate the 
system. Furthermore even if we could measure the cor- 
relations in this limit, we would not be able to determine 
h% and D^ by the procedure of Sec. vQ which relied on 
the validity of the weak-coupling approximation, which 
fails as 4> — > max . Our numerical studies only tell us 
that /if and £>f^ are increasing with packing fraction. 
In this Section we present a simple argument that in the 
limit <fi ~ ► '/'max their dependence on packing fraction is 
as 



ft? ~ 1/(1 - 0/</w) 



and that 
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1/(1 - H<h 



(29) 



(30) 



The total phase space of a finite system of hard disks 
or spheres is fractured into a number of regions ( "blocked 
states") which are mutually inaccessible. As the density 
is pushed up there are fewer and fewer blocked states 
[TT] . Eventually as the packing fraction reaches the max- 
imum for the system there is only one blocked state left. 
This can be compressed to a jammed state at </> max - The 
pressure diverges to infinity according to Eq. (fl7l). We 



random bond Ising model contains only couplings J,- o of a ,,. , ,. ,, , lU r-ssi 

.,. . , . „• , . can use this observation to deduce how hr and D h - 



^ TT 



positive sign, so it leads to a growing ferromagnetic corre- 
lation length rather than a spin glass correlation length. 
In the context of the binary mixture SOC model, if the 
ratio of the disk sizes Rab gets close enough to unity, 
we would expect that in this limit, the model would be 
in the universality class of the random field ferromag- 
net also. But this transition would be associated with 
the kind of crystal ordering visible in Fig. 13] and seems 
irrelevant to the physics of glasses. 



the effective spin Hamiltonian of Eq. Mty must vary as 
4> ~~ * 0max so as to recover the exact expression for the 
pressure of the hard sphere or hard disk gas in Eq. (17 1. 



As can be seen in Fig. 13 at low packing fraction £ 



is larger than £ and it grows with packing fraction. 
However, it does appear to saturate at around 1.5 large 
disk radii while ^sg starts to grow much more rapidly 
as the packing fraction approaches ^ ma x- This is good 
evidence that the important correlations here are spin- 
glass like and that when Rab = 1.0/1.4 the spin system 
is not behaving as a ferromagnet in a random field. The 
effective bonds generated must contain a sufficiently large 
fraction of negative bonds so that the system behaves like 
a spin glass. 

Unfortunately it is very hard to measure the correla- 
tion length £ well from simulations in the region of 
most interest, that is when <f> — > </> max . £, SG measures the 
size of the cooperatively rearranging regions and such re- 
arrangements become very slow when the required rear- 
rangements involve the cordinated motion of many disks. 
This in turn means that it takes an extremely long time 



Our argument is just a variant of the procedure of Sals- 
burg and Wood \jgj\. 

The jammed state at </> max will be isotactic to a high 
degree of approximation. That is, each disk or sphere will 
be touching z — 2d neighbors. Only a few (if any) will be 
jammed by virtue of their centers touching a plaquette 
wall and we will assume this does not occur for the state 
at 0m a x. In Sec. |XHI| a variation of the SOC model based 
on Voronoi cells is outlined where this will certainly be 
true. A finite number of disks touching the plaquette 
wall would not in any case affect the argument. Then 
in the spin mapping, such a jammed state should be a 
minimum of the Hamiltonian in Eq. Q. Suppose this 
minimum occurs at values of sf = S^". The Hamiltonian 
at its minimum, 



1 N 



(31) 



is just a constant, independent of <p with the above scal- 
ings of /if and D^ . (Here Fh" is the matrix whose 
inverse is f3D^). This expression for PT-Lmin m the par- 
tition function defined by Eq. (fn) would not then give 
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a contribution to the pressure. The pressure is actually 
determined by the contribution from the vicinity of the 
jammed state at <p max . To evaluate this contribution to 
Z let us write sf = S? + (1 — 0/</> ma x)/f ■ Because we 
are expanding about a minimum, the integrals over the 
ft are Gaussian quadratic forms in the f\~ which do not 
depend on (1 — 4> / 4> mayi ) , with our assumed scalings of 



D^- . They give a contribution to the partition function 



(1 - ^IK 



\dN 



(32) 



via the terms which comes from the changes in the in- 
tegration variable from sf to ft*. This yields Eq. (17) 



for the pressure. This result is just a consequence of the 
scaling assumed for hf and D^ with (1 — </>/0 max ). 

Note that according to this argument, both the mean 
and the standard deviation of the couplings, (which we 
generically label /3Jq and /? J respectively, without distin- 
guishing the labels p and v), will scale in the same way, 
viz as 1/(1 - </>/0max) 2 - 0h will scale as 1/(1 - </>/</> max ). 

Inserting these expressions for j3h and /3 J into Eq. (fTl) 
we recover Eqs. pi) for the correlation length. For hard 
disk systems and hard sphere systems we are therefore 
predicting that there is an actual divergence of the cor- 
relation length as <j) — > 4>max.- The circumstances where 
this behavior might be relevant to the unconstrained sys- 
tem are discussed below. 



XIII. DISCUSSION 

In a supercooled liquid, a particle is caged on time 
scales less than the alpha relaxation T a . On longer time 
scales it can diffuse anywhere in the system. In the SOC 
model, each particle is caged forever in the cell into which 
it was first inserted. 

We can measure t q from the incoherent scattering 
function: 



1 

N 



F(£,i) = -^£7,''-'""-""J 



(33) 



However, F(k,t) will never decay to zero in the SOC 
model - it will fall to a plateau and remain on the plateau 
for all time. In order to see why that happens consider 
the root-mean-squared displacement: 



rW*) = (^£R(*)-*(0)] 5 



(34) 



For the unconstrained system rMSD(t) first steadily in- 
creases with time, levels off while the particle is caged 
and finally grows to infinity. For the system under the 
SOC the cell walls ensure that r MSD {t) will saturate at 
a value determined by the size of the cell. This in turn 
ensures that F(k,t) remains non-zero for all time. 



This does not mean that the relaxation times of the 
SOC system are infinite. Consider 



c '(*) = ^E^(o)-Si(*)>> 



(35) 



and note that (sj(O).Sj(i)) = (sj) . (s*j) as t goes to in- 
finity, so in this limit C(t) approaches q. The timescale 
obtained from a study of how long C(t) takes to reach q 
would be similar to r a in the unconstrained system: the 
relaxation time r a comes about because rearrangements 
on the scale of £ in the unconstrained fluid are needed to 
relax the cages holding the particles. In the SOC model, 
rearrangements on the scale of £ are also required to al- 
low full relaxation, so the two timescales are similar. We 
leave the details to future studies. 

There is disorder present in structural glasses on the al- 
pha relaxation time scale - their molecules move so little 
that the local environment of any molecule is effectively 
disordered. However over periods of many alpha relax- 
ation times, the disorder is averaged out. Given this, 
the SOC model where quenched disorder is built in, may 
be appropriate for studying the behavior of the fluid on 
timescales of order r a . Furthermore it is from data on 
such timescales that one can obtain estimates of the cor- 
relation length £. We expect that at least when £ is large 
there is probably little difference between the point-to-set 
length scale and the dynamic length scale [I]. 

Estimates of the dynamical length scale £ are obtained 
as follows. The four-point correlation function Gi(f, t) 
defined as [48] : 

G 4 (f,t) = (p(0,0)p(0,t)p(f,0)p(f,t)) 

- (p(0,0)p(0,t))(p(r,0)p(r,t)), (36) 

should develop a plateau when the liquid starts to be- 
come glassy. The dynamic susceptibility is calculated by 
integrating G±{r,t) over volume: 



X i(t) = y I Gar.n.r'r. 



(37) 



When measured in a glassy system, X4,(i) is observed to 
grow with time, peaking at times comparable to r Q before 
decaying. As the temperature is lowered or the packing 
fraction is increased, the peak moves to longer and longer 
times (corresponding to the increase in r a ). The dynamic 
susceptibility can be thought of as a 'correlation volume' 
which reveals the scale of regions which are dynamically 
correlated [5] , providing evidence of a growing correlation 
length £ in glassy systems. 

In the binary disk SOC system, quenched disorder is 
present in the form of the random distribution of disk 
species over the cells. The growth of xsg an d £sg re- 
veal the presence of growing amorphous order. Because 
for the SOC system the quenched disorder persists for all 
time, not just for timescales up to r Q , if Xi(t) were mea- 
sured in the SOC model it would grow and then saturate 
at xsg- 
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It is our belief that SOC models of hard disks and 
spheres can therefore describe the increase of £ with pack- 
ing fraction, at least as regards the value of the exponent 
v. We do not expect the value of ^> max to necessarily coin- 
cide with the packing fraction of the divergence in the un- 
constrained system - after all, </> max would be of slightly 
different value if the cells had not had a square shape 
or even were of random shape. A way of constructing 
"random" cells would be to equilibrate the binary disk 
or sphere system and then use as the cells the Voronoi 
cells of a single configuration as the cells. Because of 
this built-in randomness, this same procedure could be 
used to model the striking glassy features of monodis- 
perse spheres. (For the Vor onoi cell version of the SOC 
model, the argument in Sec. 
the other hand, for such eel 



XII is clearly exact 30J . On 



s it would be impossible to 



IVI. 



carry out the analytical calculations in Sec. 

The divergence of £ as 4> ~ * </>max is likely to be accom- 
panied by a divergence of the relaxation time of the SOC 
model. Note that such behavior is not that expected of a 
G point [49]. At a G point £ and r both diverge, but the 
pressure remains finite. At </> max the pressure is infinite, 
as it is also a jammed state. 

Our value for ^ max is quite close to the estimates of the 
value of the random close-packing fraction: (f> rcp m 0.84 
|33j . We think that this similarity is not an accident. 
Both the packing fractions, <f> max and rc p, are obtained 
from situations where the phase space of the hard disks 
has been curtailed so that the system cannot stray far 
from its initial state. </> max will depend on the choices 
made for the shape of the cells. It will also depend on 
how the large and small particles are assigned to the cells. 
In our work this has been done randomly but one could 
build into the distribution if desired the local correlations 
of the unconstrained system. It is also known that the 
random close packing fraction <p rcp is not well-defined: it 
has a small dependence on the protocol used to determine 
it [33J. 

When studying the unconstrained hard sphere or hard 
disk system, some protocol has to be adopted to see 
glassy behavior, such as a finite compression rate, and 
this will result in the pressure going to infinity at some 



packing fraction less than that of the densest crystalline 
state. In true equilibrium, the pressure of course remains 
finite unless the system is at the maximum density of 
the crystalline state. We believe that the glass state is 
well-defined provided that the alpha relaxation time r a 
is such that 1/t q is greater than the rate for phase sepa- 
ration and crystallization in the case of binary mixtures, 
or the time scale for crystal nucleation and growth gen- 
erally. A finite compression rate should not modify the 
quasi-equilibrium approach to glasses (like that in this 
paper) provided that it is small compared to 1/t q . Since 
the alpha relaxation time r Q is expected to grow with £ 
as In t q i-v £v [SD] , then for a fixed compression rate one 
can only hope to obtain the growth of £ up to a com- 
pression rate determined value. But within these various 
constraints we believe the glass problem is well defined 
and that SOC models are a useful way of studying some 
aspects of it. 

We would expect the SOC model of hard disks or 
spheres to be most useful at densities above (f>d, the 
density at which timescales increase rapidly. This den- 
sity can be quite well-understood with the aid of mode- 
coupling theory. For hard spheres <pd « 0.58 and for hard 
disks 4>d ~ 0.78 11 . In the case of hard disks in the SOC 
model, timescales were seen to increase very rapidly at 
a rather similar density. This is because at such den- 
sities the timescales are long because they involve col- 
lective rearrangements of the disks on a length scale £ 
and collisions with the walls of the plaquette are becom- 
ing insignificant. Alas, this very rapid increase makes 
numerical investigations at densities above (f>d very chal- 
lenging. 
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